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Abstract 

We describe some detailed numerical simulations of Disoriented Chiral Condensates 
(DCCs), using the chiral lagrangian as a controlled long- wavelength description. We 
focus on the possibility of multiple, independently coherent domains, and investigate 
the degree to which the DCC signal is attenuated. As an intermediate step in our 
analysis we compute the expected number of detector events in each isospin and mo- 
mentum channel for a given asymptotic classical field configuration. We find that for 
sufficiently large initial field strengths, the non-linear interactions between domains 
become important and can lead to a randomization of isospin orientations. Neverthe- 
less, we argue that viable signals exist for DCC detection, even in the case of multiple 
domains and strong domain-domain interactions. We briefly discuss some long-lived 
'pseudo-bound state' configurations which arise at large field strengths and might be 
observable in HBT correlations. 
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1 Introduction 



One of the most exciting possibilities at RHIC is that collisions of heavy nuclei might lead 
to coherent, pion field configurations which are produced by the disorientation of the chiral 
condensate (DCC) [l], ||, 0, ||, ^ (for recent review articles, see || 0). Such coherent con- 
figurations could lead to large, characteristic fluctuations in the number and type of pions 
measured in detectors. Unfortunately, the evolution of quantum fields in the aftermath of 
the heavy ion collision is a difficult problem, involving both strong coupling effects and non- 
equilibrium statistical mechanics, and remains to be solved despite much theoretical effort. 
Most previous work has been performed in the context of the linear sigma model, sometimes 
using an approximation like the 1/N expansion to control the quantum corrections p|-|l9|. 
However, the relationship between linear sigma model dynamics and real QCD dynamics is 
unclear. 

Our approach here will be somewhat different. We adopt an agnostic approach to the 
formation of DCCs, and simply try to characterize their evolution at subsequent times, with 
an eye toward experimental signatures. We make a detailed investigation of the phenomenol- 
ogy of DCCs using the classical field equations derived from the chiral lagrangian (non-linear 
sigma model) to describe their evolution from initial domain (s) to the final state which arrives 
in the detector. The field equations of the chiral lagrangian can be organized in a momen- 
tum expansion, and therefore yield a controlled approximation for sufficiently slowly-varying 
pionic configurations. While this approach has no advantages for the formation problem, 
it does provide a way of simulating, in a controlled fashion, the subsequent evolution of 
domains once formed. Indeed, quantum corrections (loop effects) are suppressed as long as 
the configurations simulated are sufficiently "soft" relative to the scale &7rF n — 1 GeV. We 
verify that this treatment is self-consistent: initially soft configurations do not lead to the 
build up of hard sub-configurations which cannot be treated within the chiral lagrangian. 

We will be particularly interested in the possibility of multiple coherent domains, since it 



seems unlikely that a single large domain would remain in the wake of the collision [13], [14 



17]. Therefore, we simulate the interactions between domains with different initial isospin 
orientations. For sufficiently weak fields, interactions can be neglected, leading to particle 
distributions which are simply linear superpositions of those expected from each separate 
domain. In this case the effect of iV multiple domains is mainly "statistical", leading to a 
narrower width (proportional to 1/yN) in fluctuations. However, at sufficiently large field 
strengths, which we quantify below (but which are well within the roughly expected range), 
the interactions are significant and have a strong effect on the field configurations, even if the 
number of domains is small. Our results are discussed in detail below, however the general 
observation one can make is that non-linearities can degrade the DCC signal, leading to a 
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somewhat narrower distribution in fluctuations of the charge ratios. Despite this, we find 
that promising signals of multiple DCC domain formation remain even in the non-linear 
cases. We also briefly discuss some interesting 'pseudo-bound state' (PBS) configurations 
which appear at large field strengths. These configurations have particularly long lifetimes 
and might be observable in HBT correlations at RHIC. 

The organization of this paper is as follows. In section 2 we review the chiral lagrangian 
approximation to low-energy QCD, and derive the corresponding classical field equations. In 
section 3 we discuss various aspects of multiple domains, including the straightforward sta- 
tistical analysis. In section 4 we describe our simulation and our numerical results, including 
details of the non-linear behavior. In section 5 we present our conclusions. Appendices A, 
B and C contain some useful formulae, as well as a description of our "detector" subroutine 
and of how isospin rotations of our configurations were performed. In appendix D we give 
additional details concerning the simulations and the parameters used. 



2 Chiral Lagrangian Description 

The chiral lagrangian is the most general effective lagrangian consistent with the symmetries 
of low energy QCD [pOfl . As such, it accurately describes the soft dynamics of pseudo- 
Goldstone mesons, including the virtual effects of heavier particles such as baryons, the rho 
meson, etc. which have already been integrated out. We can classify the possible terms by 
the number of derivatives, with the first term given by 

C = ^-trid^d^E) (1) 

where 

S = e w/i ^ . (2) 

The effects of higher order terms are proportional to powers of the momentum, as are any 
quantum loop effects calculated with ([]]). This implies that for sufficiently soft configurations, 
the classical field equations derived from (jl]) describe the full quantum mechanical evolution. 
For this reason, we believe that the chiral lagrangian is superior to the usual linear sigma 
model for studying the dynamics of DCC's, at least after formation. Our approximation 
relies on a momentum expansion, rather than 1/N or Hartree-Fock. 

We define a four component real field <fi = (a, <fi) such that F % H = ol + ir ■ with the 
condition 2 = F% everywhere. It is clear that we can do this, since 

e ir-A = cogfl^i)/ + % ZlA sin(|A|) . (3) 
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In terms of 0, we have 

C = l -d^ ■ , (4) 

with the condition 2 = F 2 . 

Rather than deal with the complicated constrained equations which result from (f|), we 
relax the condition on <fi 2 slightly by replacing it with a potential term in the lagrangian 
whose minimum falls at 4> 2 = F 2 . This yields the linear sigma model with potential 

V{4>) = \[4?-F%\ (5) 



and lagrangian 



C = l -d^ ■ - V(<j>) = ■ - - F 2 } 2 . (6) 



The limiting case A —>■ oo yields the original non-linear sigma model. 

In reality, the chiral symmetry is slightly broken by bare quark masses, resulting in 
non-zero pion masses. Taking this into account yields 

£ = ■ M - ~ ~ ^ ~ + Ft). (7) 

This lagrangian yields the classical equations of motion 

d 2 4> = -{\(<p 2 -F 2 )+ml]J (8) 
d 2 a = -[X(<p 2 ~F 2 )+ml]a + m 2 n F n . (9) 

In the non-linear sigma model limit we wish to consider, a is just F n plus a very small 
quantity. Because floating point variables do not deal well with this arrangement, it is more 
useful to redefine a — > a + F n , with equations of motion 

d 2 $ = -[X(4> 2 + 2F w a) + ml\4> (10) 
3 2 ~o = -[\(4> 2 + 2F n d)+ml]a-\{<t> 2 + 2F 7T a)F 7T . 



3 Multiple Domains 

A telltale sign of DCC formation would be large fluctuations in pion isospin abundances. 
Under incoherent production of hundreds or possibly thousands of pions, one would expect 
nearly equal numbers of each type to be emitted, whereas with coherent DCC domains, the 
abundances would be a function of the domain isospin orientation. If one domain happened 
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to point in the ttq direction, then each of the outgoing pions would be a 7To, an improbable 
occurrence under incoherent production. In the case of a single, coherent domain, the exact 
probability distribution for 

/ s „ . . „ (id 



is easily seen to be dP(f) = df/2y/J . 



The case of k multiple domains is more complicated, but in the low field strength ap- 
proximation it can be simplified by the neglect of interactions, so that the distribution of 
pions is the sum of the distribution from each single domain. If the size of each domain is 
taken to be roughly the same, and hence total particle output iV = N n o + N w + + N n - does 
not differ from domain to domain, then the total / is simply the average of that for the 
individual domains, 



total ^— , „ T TAT 7 / j J 



(12) 



£ Nj kN k 

As was also shown in f2~I|| , the distribution for this total / can now be computed as a func- 
tion of k. Analytical results are difficult to obtain, but by randomly generating domain 
orientations (with a flat probability distribution on the three sphere), the probability dis- 
tribution curves can be calculated numerically. By generating large quantities of samples, 
the error bars in the distribution curves were lowered enough to be negligible in the graphs, 
approximating the limiting case of a continuum probability distribution. Figure [l| shows the 
behavior of P(f) for increasing values of k. 



ir 1-5 




Figure 1: Probability distribution of / [P(/)]for varying k 
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The general behavior is a result of the central limit theorem, which tells us that as 
statistical sampling is performed repeatedly, the average value of the resultant variable stays 
the same, while its standard deviation is proportional to l/y/k. The proportionality constant 
is easy to fix by calculating the standard deviation for the k = 1 case, whose probability 
distribution was given above as dP(f) = df/2y/J. This yields 



^coherent = V / (/ 2 ) " UV = ^ • (13) 



We can repeat this calculation for the incoherent case with N particles emitted, using 
the probability distribution for single particle production, 

p(f) = ls(f) + \s(f-i)- 



In this case 



^incoherent 3 V iV ' \~' 



While we do not know precisely how many domains (if any!) are likely to be formed in a 
typical collision, it is difficult to imagine that k is much larger than of order 10 or so, simply 
due to the restricted geometry of the collision region. Therefore, N is likely to be large 
compared to k and there is still a significant qualitative difference between the distribution 
of results from incoherent and coherent production. The overall level of fluctuations in / 
should be larger in the coherent case. Whether this difference could be observed at RHIC 
depends of course on the probability of DCC formation per collision. 

The coherent analysis would be modified if some fraction a of the outgoing pions were 
assumed to be produced by incoherent background processes. For (N^ , N^ + , _) coher- 
ently produced and (N" Q , N" + , N"_) incoherently produced pions, (where N" otal = a(N' total + 
N total)) i the value of / is (neglecting effects of order 1/y/N") 

K + K 1 

f = ~W' ° _j_ J\77/ = (1 _ a ) /non-corrupted + 77^ • (15) 
JV total JV total 6 

Using this the standard deviation of / becomes 

Ocohcrcnt = (1 - a) ^— - . (16) 

3v5/c 



Hence incoherent corruption narrows the spread in / around 1/3. 

In the above discussion we ignored the interactions between the coherent configurations. 
Clearly this is justified for sufficiently weak fields, but when the pion field strengths are 
large the interactions between different domains can be significant. It is well known that 
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the solution of classical field equations is equivalent to the resummation of all tree graphs in 
perturbation theory. Nonlinearities in the classical field evolution therefore represent large 
rescattering effects due to tree-level vertices. As mentioned previously, all loop graphs in 
chiral perturbation theory contain extra powers of the external momentum squared, and are 
hence suppressed. We expect rescattering to be important when the pion number densities 
are sufficiently large. The effect of the interactions is not clear a priori, although one might 
argue that they tend to randomize the isospin orientations and lead to smaller, not larger, 
fluctuations in / for a fixed number of domains. To gain some insight on this question, 
we performed numerical simulations of domain-domain interactions using the low energy 
effective field theory described in the previous section. 

In our simulations we allowed the field strengths to vary in the range of \tt\ ~ F^. Due to 
our lack of control over the formation process, we are unable to compute with any precision 
the expected number or energy density of pions in the DCC. However, a simple estimate 
based on a critical temperature of chiral symmetry restoration T c ~ 200 MeV, and the 
assumption that the energy density of pions in the DCC is roughly that of the hadron gas at 
T c , yields field strengths as large as (4 — 5) F n . As we will see below, the evolution of multiple 
domains is significantly non-linear at field strengths which are as small as \tt\ ~ lAF n . 

4 Simulation and Results 

We performed simulations of domain-domain interactions in a cylindrical geometry with ax- 
ial symmetry. The domains themselves are cylinders placed end to end along the beam axis, 
with domain walls in between. This is of course an idealization of the actual configuration, 
chosen for simplicity and to make our simulations tractable. 

4.1 The Difference Equation 

A computer simulation written in C++ was used to evolve eq. (0). Imposing axial symme- 
try, the left side of the equation becomes 



d 2 cP 



d 2 <p d 2 <p 19 d(f>, d 2 <f) d 2 <p 




1 d<p d 2 <p 



(17) 



dt 2 dz 2 r dr dr dt 2 dz 2 



Defining a discrete function 





and using the following finite-difference approximations 

dq ~ 2Aq 



(19) 
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d 2 <f) 



<f) i+ i - 2<j)j + <f)j_i 
Aq 2 



(20) 



dq 2 



the equations of motion 



become 



'i+l,j,k — 




(21) 



(Here Ar = Az, and G^j^) is the expression on the right hand side of the original equa- 
tion). 

By storing the two previous time slices of the field {4>i,j,k and the newer value 

(<f>i+ij,k) can De calculated, allowing for the time evolution of the fields given an initial 
configuration and its time derivative. 

4.2 Initial Configurations 

An initial configuration is described by its number of domains k, their positions, sizes, field 
strengths, and isospin orientations. For simplicity we assume the size and strength of each 
domain is the same (owing to a similarity in formation processes). We also choose the field 
configurations of the first two timeslices to be the same. As a final simplifying assumption we 
have ignored the initial expansion velocities of the domains. Such velocities would probably 
decrease the nonlinear effects, since the pions would disperse more rapidly to low densities. 

Distinct fixed values of k are treated differently. 



Here the initial conditions for each run are characterized by one vector v representing 
the orientation and strength of the single domain in isospace. 

Varying the field strength in the k=l case yields no change in the distribution curve for 
/. It can be easily shown (e.g., using the results in appendix C) that the final values of N no , 
N n+ , and N n _ are proportional to the square of the three components of v: 



The proportionality constant cancels out in the definition of /, leaving it and its statistics 
invariant to the initial strength. 



The initial configurations for the k=2 case can be parameterized by two isospin vectors 
vi and V2- To generate the P(f) curves, the field strength of the two domains was fixed at 



k=l 




(22) 



k=2 
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the given value \vi\ = \v%\ = s, multiple runs with random orientations were performed, and 
the resulting /'s were binned, much like in the weakly interacting case. 

Each run requires tens of minutes, hence it is too time consuming to perform the O(10 9 ) 
simulations needed to smooth the probability curves. Instead, the rotational techniques 
developed in appendix C were exploited, and the simulation was only run for v\ and v-i 
orientations in the 7r — 7r + plane, symmetric around the 7r axis. Rotating this representative 
class yields all possible initial configurations. 

The actual functional shape of the initial field configurations is specified by a magnitude 
and angular dependence. The magnitude is given by 

\n(r,z)\ = s*g{\z\)g(r) , (23) 

with g(x) = 1/ (l + e kl (' x ~ a ^ . This function (for appropriate kis) is nearly constant for 
x < a, drops of suddenly at x ~ a, then stays approximately for x > a. The angular 
dependence, measured in the ttq — tt + plane from the 7To axis, is 

2 a 

ln(r,z) = arctan(/c 2 z), (24) 

7T 

which plateaus to +a for z > and —a for z < 0. These equations describe a cylinder 
centered at z = r = subdivided into az>0 part with isospin in the +a direction and a 
z < part with isospin in the —a direction. Values of k\, &2, an d a were chosen to fix the 
physical dimensions. Each isospin region has a height and diameter of 10 fm, with a 'skin' 
thickness of 1 fm. (The 'skin' is defined as the region where the field amplitude drops from 
90% to 10% of its maximum value.) Additional details of the computational techniques used 
in the simulations can be found in appendix D. 

Figure |2| gives sample time snapshots in the evolution of the pion fields, where light white 
(dark black) corresponds to the maximum (minimum) value on the field in the given snap- 
shot. Runs were continued until the Klein-Gordon Number Operator converged, signaling 
that the field had spread out far enough that the nonlinear terms had become negligible. 

The progression of P(f) for various initial field strengths was constructed and is shown 
in figure |^. The dotted line shows the low field strength approximation. Successive curves 
are then shown for initial field strengths of O.OIF^., 0.5F n , 1.0F n , l.3F n , lAF n , 1.5F W , and 
l.dF^. The trend towards a more randomized probability distribution is apparent, although 
the width of the distribution and the area under the tails (in particular, / ~ 1) is not 
changed significantly. One way to understand this is to note that in order to get an event 
with / close to unity, the initial domains must all be oriented in the same direction. This 
is of course less and less probable as the number of domains k gets larger, but once the 
domains are formed in this way the common orientation suppresses their interactions - in 
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time = 0.0 fm time = 5.7 fm time = 16.9 fm time = 29.8 fm time = 41 .9 fm 
Figure 2: Numerical evolution of the fields 

other words they act like one big domain. Thus we expect the number of / ~ 1 events to 
be preserved even at large field strengths. 

k > 2 

Computational power limited our ability do do an exhaustive analysis of the k>2 case. 
For the k=3 case, initial configurations are parameterized by the three unit vectors and an 
overall strength. Even exploiting the rotational symmetry and fixing field strength, three free 
parameters remain. Building a representative class of configurations, sampling each variable 
at S values would would need at least S* 3 computer runs. Even the marginally acceptable 
value of S — 10 would require 1000 runs, which at greater than 2 hours per run on our DEC 
Alpha Workstation would require almost 3 months. Although more powerful computers 
could be assigned to this task, we feel that the important qualitative features of the domain 
to domain interaction are already featured in the k=2 case. This intuition is based in part on 
the (explicitly confirmed) observation that when two domains are separated from each other 
at larger distances, interactions become quickly negligible. Thus, in the case of k domains 
arranged colinearly (as required by the axial symmetry) only the nearest neighbor domain 
walls would interact significantly. Of course this single domain wall was well studied in the 
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Figure 3: Probability distribution of / [P(f)) for varying field strength 

k=2 case. 

4.3 Nonlinearities and Pseudo-Bound States 

Increasing the strength of the field configurations qualitatively changes the evolution of 
the fields. We carefully examined the runs to determine the critical amplitude at which 
nonlinearities set in. In these observations a surprising behavior was noticed. For values of 
initial strengths between and 1.2.FV, evolution proceeded similarly to the non-interacting 
case. However, for slightly higher field strengths we observed that a large fraction of the 
energy bunched up in the central region and held together for long periods of time. Figure |] 
shows the energy in a 10 fm spherical region centered around the origin vs. time. As can be 
seen in these graphs, once the threshhold has been crossed a pseudo-bound state forms. 

This phenomenon was found to be robust: it was replicated in runs in which the initial 
conditions were varied. We also studied the evolution of spherically symmetric configura- 
tions using separately written code. Again, once a certain initial strength was crossed, the 
configuration held together in a pseudo-bound state for much longer periods of time. It is 
in fact the long evolution times required for the configuration to completely disperse that 
limited our ability to compute the probability distribution in / beyond the 1.6.FV case. 

These much longer dispersion times are potentially observable in HBT (Hanbury Brown- 
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Figure 4: Central energy vs. time 



Twiss) correlations^. (For a recent review of these techniques, see pffi| .) Pseudo-bound state 
dispersion times as long as 10~ 22 s (versus 10 _24 s for the noninteracting case) were observed, 
and we have every indication that longer lifetimes are possible. We intend to report in more 
detail on PBS dynamics in future work. 

We note that although the configurations evolved very differently for higher field strengths, 
the P(f) curves were only slightly modified. This fact leads us to assert that, at least up to 
the strengths that we were able to study, the noninteracting case remained a good approxi- 
mation for the full evolution. 



5 Signals 

In this section we discuss two signals which can be used to detect or exclude DCC production. 
We focus on simple hypotheses with specific probabilities pk for producing k domains of some 
fixed characteristic size. 



We thank K. Rajagopal and D. Rischke for pointing this out to us. 



12 



5.1 Width Measurements 



Purely incoherent pion production should yield a spread in values of / equal to a = §y j| 
(equation |1^). Any significant excess of this value would be strong evidence of DCC forma- 
tion. Although the range of probable RHIC pion yields is presently unknown, experience 
from previous heavy ion collisions would suggest values in the hundreds or thousands, yield- 
ing a = 0.01 to 0.05, a range that would be surpassed for up to k=35 domains (equation |T3|). 

In reality, k domain DCCs would be expected to form in only some fraction pj, of the 
events, yielding a smaller value of a. Data can place upper bounds on each value of pk from 
the following analysis. Assuming that only k domain DCCs and incoherent events form (for 
some given k), the probability distribution curve for / becomes 

P{f) = Pk -Pk- domains^) + (1 ~ Pk) -^incoherent (/) > (25) 

hence 

""measured 

/ df(f~ < f >) 2 [ Pk P k 

—domains 

(f) + (l- Pk ) P mcohcrent (/)] (26) 



Pk "coherent + i 1 ~ Pk) O? 



incoherent 



Using equations O and O, and reintroducing incoherent background corruption discussed 



in section 0, we can obtain an upper limit of pk 



5 

Pk < -k 



9Na 2 m 



2(1 -a) 2 N -5k 



(27) 



Figure |5| shows a plot of numerically calculated values of width vs. initial field strength 
(for the k = 2 case). Despite the level of nonlinear activity that occurs at higher strengths, 
the width does not change more than 10%. Up to the range of simulated cases, the noninter- 
acting case is a good approximation. For higher strengths, there is probably a trend toward 
somewhat smaller widths. 



5.2 Counting Rare Events 

Another feature of incoherent production is that the distribution of / peaks sharply at 1/3, 
and is nearly zero elsewhere, yielding an exponentially small tail of events where / is close 
to either or 1. From elementary statistics the distribution of / is 
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Figure 5: Standard deviation vs. initial field strength 



P(f = i/N) 




(28) 



where i is restricted to being an integer from to N. Setting N=100 and integrating over 
the range 0.6 < / < 1 yields the small value of 4.3 x 10 -8 . Thus, given (for instance) 10 7 
events, there is still only a 43% chance that / > .6 be measured even once. 

A similar anaylsis for the coherent case was performed. For instance, with one domain, 
P(f) = —7=, for which 23% of the events lie above / = 0.6 (a much larger 2.3 x 10 6 of 10 7 

assumed events). Now, assuming that only the incoherent and one domain coherent cases 
contribute, the fraction of rare events t would be 



t=(l- Pl ) t' + Pl t" 



(29) 



where if is the number of rare events expected from the incoherent case, t" is the number 
from the coherent case, and p\ is an element in the set of p^s introduced in the previous 
subsection. Presumably more than just these two cases would exist, so this equation sets an 
upper limit on p\ 



Pi 



< 



t-f 
t" - V 



(30) 



In this manner, using the probability curves of / calculated in previous sections, data for 
other values of k were calculated. Table [l] shows the expected fraction of tail events for given 
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Number of domains k 


Fraction of events 
with / > 0.6 


Upper bound on p k (< jnSjj) 


k — 1, all field strengths 


2.3 x 1CT 1 


4.3 £- 1.9 x 10~ 7 


k — 1, all field strengths ' 


1.8 x 1CT 1 


5.4 £- 2.3 x 10~ 7 


k = 2, Low field strengths 


1.2 x 1CT 1 


8.3 £ - 3.6 x 10~ 7 


k = 2, Low field strengths t 


7.4 x irr 2 


14 t - 5.8 x 10" 7 


k = 2, High field strengths 


1.1 x 10" 1 


9.1 t - 3.9 x 10~ 7 


k = 2, High field strengths t 


7.1 x icr 2 


14.1 t — 6.1 x 10" 7 


fc = 4 


5.0 x 10~ 2 


20 £-8.6 x 10" 7 


k = 4 t 


1.8 x 10" 2 


56 £-2.4 x 10" 6 


fc = 8 


8.4 x 10~ 3 


120 t- 5.1 x 10~ 6 


fc = 8 f 


1.4 x 10" 3 


710 £-3.1 x lO" 5 


fc = 16 


3.4 x 10~ 4 


2900 £ - 1.3 x lO" 4 


fc = 16 f 


1.1 x 10" 5 


9.1 x 10 4 £ - 3.9 x 10" 3 


fc = 32 


7.7 x 10" 7 


1.3 x 10 6 £ - 5.9 x 10- 2 


fc = 32 t 


8.8 x 10~ 10 


2.4 x 10 7 £ - 1.0 



(t = with 20% incoherent background) 



Table 1: Rare events 

domain configurations (varying k, fraction of incoherent corruption, and field strengths for 
k=2). Also included are upper bound expressions for the p^s. 

From this data, we can conclude that the fraction of events lying in the tail falls off 
exponentially with fc. Figure |6| shows a fit to the form 

P t ^ = Ae- Bk . (31) 

with fitted values A = 0.147 and B = 0.370. From this, it is estimated that k could increase 
to 40 before the number of tail events would be comparable to the incoherent case. Similarly, 
with 20% incoherent background included, A = 0.20, B = 0.61, and k can still increase to 
around 25 before the number of tail events would be comparable to the purely incoherent 
case. 

6 Conclusion 

We have conducted a detailed numerical investigation of multiple domains of disoriented 
chiral condensate. Although the dynamics of domain formation in the wake of a heavy ion 
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Figure 6: Tail events versus k 

collision is far from completely understood, it is plausible that multiple domains are more 
likely to occur than a single coherent domain. We have concentrated on the evolution of 
multiple domains after formation, using the chiral lagrangian as a long-wavelength approxi- 
mation for the dynamics. Our simulations follow the pion fields out to asymptotic distances, 
where we can characterize the number of particles per isospin and momentum bin that would 
be detected in an experiment. In choosing the initial conditions for our configurations, we 
opted for simplicity in neglecting important factors such as the expansion velocities. Given a 
particular dynamical model of the formation process, this could easily be incorporated into 
our simulations. 

For small amplitudes, corresponding to low number densities of pions in the semi-classical 
fields, a simple statistical analysis suffices to characterize the isospin fluctuations of multiple 
domains (see figure 1). However, as the field strengths are increased, we see ample evidence 
of nonlinear behavior and significant interaction between domains. These nonlinear effects 
are already present for field strengths of order ~ F w . For the case of two domains (k — 2), 
we were able to compute the effects of interactions on the probability distribution of isospin 
fluctuations (figure 2). While the shape of the curve is changed somewhat from the weak-field 
case, the probability of rare events (in particular, where f is close to unity) and the width 
of the distribution are not strongly affected. We suspect that the rare events with large f 
will prove to be the most robust signal of DCCs, surviving even multiple domain formation 
and large field strengths. Some interesting phenomena were observed in the evolution with 
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high initial field strengths. Relatively long-lived 'breather' or 'pseudo-bound state' config- 
urations were found, in which the pions remained trapped on timescales many times larger 
than the light-crossing time of the initial configurations. In future work we hope to under- 
stand their dynamics better, as well as the prospects for their detection through correlation 
measurements. 

Computational limitations prevented us from performing complete investigations of the 
case of k > 2. While we can perform individual runs, the time required to compute the 
associated probability distribution for / was prohibitive. Nevertheless, for the simplified 
cylindrical configurations we considered it is plausible that the k = 2 case captures the 
correct qualitative behavior for larger numbers of domains, as in any individual run the 
evolution is seen to be dominated by nearest domain interactions. Our results suggest that 
even events with a large number of small domains are capable of providing a detectable 
signature. We are hopeful that analysis of future RHIC data will allow stringent upper 
bounds to be placed on the probabilities of multiple domain formation, or alternatively the 
discovery of coherent phenomena in heavy ion collisions. 

The authors would like to thank John Harris, Rudy Hwa, Vincent Moncrief, Krishna 
Rajagopal and Dirk Rischke for useful discussions and comments. Nick Evans and Stephen 
Selipsky are acknowledged for early participation in this project. This work was supported 
in part under DOE contracts DE-AC02-ERU3075 and DE-FG06-85ER40224. 



7 Appendix A : Relations Between n and <f> Fields 

The and n fields are related via 



(32) 



Using the identity 



,iA-f 



cos(|A|) + i 



(33) 



we can write out the above equation 




(34) 



which implies 




TT 



(35) 
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and 



F^sin(^) . (36) 



By inverting the equation for 0o above, we get 

\ir\ = F n arccos(-^-) , (37) 

and since tt and <fi point in the same direction (as seen from the equation solving for above), 
we get 

tt = F 7r arccos(— ) — - . (38) 
F^ |0| 

However, this will not suffice for the case where the values of the fields are approximated. 
The problem is that there are four fields while there are only three tt 's. The reason that 
this is possible is because there is an implicit condition on the 's - that 2 = 0q + <p 2 = 1. 
If this is broken even slightly, it is possible that O > F w , and then cos _1 (|r-) becomes 
imaginary. 

To avoid this problem, a fourth parameter is coupled with the 7r's, a K field that measures 
the length of (f) (K 2 = <p 2 ). Using this, equations ^ and [SB] become 

0o = ^cos(^) (39) 

= i^sin(i^) . (40) 



7T 



Inverting these yields 



and 



K = \<S>\ , (41) 



tt = F n arccos(— )— . (42) 



Now we can assume that K ~ F n and throw away the information in K. 



8 Appendix B : The Particle Detector 

The number of particles can be detected by using the usual Klein-Gordon number operator 
in the limit of the simulation where the wave is dispersed and the nonlinearities have become 
negligible. The number operator is 

N = j^t- 3 J d 3 p 4 op . (43) 
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As we are in the classical limit, we can treat all these operators as c-numbers and simply 
solve the above in terms of the fields. To do this we start with 



tt(s) = / e* >* + a\ e* -) . (44) 



After some algebra, this equation can be inverted to solve for a, yielding 



V2 



d x 



ir(x) + i 




TtlX) 



-l p-x 



(45) 



V2 



e iE ? FT r - 



Ep tt(x) + i 




71" UE 



(46) 



It might seem natural to simply insert this into the definition of N and derive an explicit 
expression. However complexities arise due to the non-local nature of the N operator. Hence, 
we separately calculate the ap field then use this to obtain N. 

Two useful equations are those of the Fourier transform in axial and spherical coordinates, 
given by the spherically symmetric case 



1 r°° 

FT(P r )[f(r)) = I d 3 x f{r)e~ l x ' p = 4vr— / dr f(r) r sin(rP r 

P y J 



(47) 



and the axially symmetric case 
FT(P r , P z )[f(r, z)\ = I d 3 x f(r, z) e~ l = 2vr 



oo /-oo 



dz dr f(r) r J (rP r ) e~ l zP * . (48) 



J-oo 



In addition to its use in calculating the number of particles, transforming to momentum 
space allows us to inspect against possible aliasing problems. Because of the axial (or spher- 
ical) symmetry and complicated shapes involved, this would be difficult to do analytically. 
However, a plot of a p vs. \p\ (for representative values of p) can be visually checked. As in one 
dimensional Fourier Transformations, discretization copies the primary Fourier transformed 
image centered around zero momentum in a repeating pattern at higher frequencies (albeit 
deformed). By viewing a p we insure that each copy decays fast enough as to not interfere 
with its neighbor. Figure shows a sample a p , which can visually be seen to decay before 
the second image (which is not shown). 
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Figure 7: An example of typical values for a p vs. \p\ 

9 Appendix C : Rotating the Configuration 

To minimize the number of simulation runs, we used the results of one run to calculate the 
detected pion distributions resulting from initial conditions which differ from those of the 
first run by a rotation in isospin space. The definition of the Klein-Gordon number operator 
is 



N =T^l d *P4 a v ■ (49) 



In our case we have three associated number operators, for 7t 0i+i _: 



' = (2^/^ ( < :)t < : • (50) 

An isopsin rotation on the configuration yields a new set of a/s via the usual transfor- 
mation, 

a j = R), a?' , (51) 

where R is a 3 x 3 5*0(3) matrix. Applying this to the iVj's, one obtains an expression 
involving cross terms (a^a k , so a new tensor is introduced 

N*> k = j^y 3 J <?p (4)t a} , (52) 
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where the single indexed N l, s are just the diagonal elements of the double indexed N^ ,k , 
Now the N's can be trivially rotated via the usual tensor rotation law, 

W' k = R), R\, N^ k ' . (53) 



10 Appendix D : Computational Technique 

Besides initial configurations, each numerical simulation had a series of internal parameters, 
whose values were independent of the physics of the problem but needed to be determined 
in order to insure proper results were obtained. These parameters included the spatial 
discretization Ax, time step At, total evolution time T, and linear sigma potential scaling 
A. 

Of these, Ax and At proved to be quite simple to fix, as the final results were insensitive 
to the chosen values, aside from cases where evolution became obviously unstable. Chosen 
(Ax, Ai)'s within the stable region yielded final values which agreed within reasonable er- 
rors, so for a given Ax, the coarsest stable value of At was chosen, thus decreasing needed 
computational power. It should be noted that the required At decreased as field strength 
and A were increased to high values, because when this is done, higher frequency oscillations 
become more dominant in the evolution. 

A reasonable value of A was found by holding everything else fixed, then calculating and 
plotting the number operator for varying A (see figure §). From this plot we were able to 
visually pick out the convergence point of the curve. The curve converges for values of A 
in the hundreds, but as the computational cost of raising A turned out to be low, the much 
more conservative value of A = 6000 was chosen. 

The remaining parameter, evolution time T, proved to be the least straghtforward to 
work with. In principle the same technique was used as for A, but in practice difficulties 
resulted due to the large and sudden variation in decay time for differing initial conditions 



(see section [O]for details of this phenomenon). Longer evolution times stretched the limits 
of our machines' memory and processor speed, and in the end were the deciding factor in 
limiting initial field strengths less than or equal to 1.6 F n . Figure || shows a typical plot of 
calculated N versus T illustrating how the convergence process took place. 
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